% Solve the model steady state 
% .. `bar' denotes steady state values

ibar = 1/betta*pibar;   % Nominal interest rate
i_annbar = 4*ibar;      % Annualized nominal interest rate
pi_annbar = 4*pibar;    % Annualized inflation rate
abar = 1;               % Technology

% Real marginal cost
Omega_21 = ( ((1-theta_2*pibar^(eta-1))/(1-theta_2)) / ((1-theta_1*pibar^(eta-1))/(1-theta_1)) )^(1/(eta-1));
Psi_21   = ((1-theta_2*betta*pibar^(eta-1))/(1-theta_2*betta*pibar^(eta))) / ((1-theta_1*betta*pibar^(eta-1))/(1-theta_1*betta*pibar^(eta)));
Omega_31 = ( ((1-theta_3*pibar^(eta-1))/(1-theta_3)) / ((1-theta_1*pibar^(eta-1))/(1-theta_1)) )^(1/(eta-1));
Psi_31   = ((1-theta_3*betta*pibar^(eta-1))/(1-theta_3*betta*pibar^(eta))) / ((1-theta_1*betta*pibar^(eta-1))/(1-theta_1*betta*pibar^(eta)));
Omega_41 = ( ((1-theta_4*pibar^(eta-1))/(1-theta_4)) / ((1-theta_1*pibar^(eta-1))/(1-theta_1)) )^(1/(eta-1));
Psi_41   = ((1-theta_4*betta*pibar^(eta-1))/(1-theta_4*betta*pibar^(eta))) / ((1-theta_1*betta*pibar^(eta-1))/(1-theta_1*betta*pibar^(eta)));
Omega_51 = ( ((1-theta_5*pibar^(eta-1))/(1-theta_5)) / ((1-theta_1*pibar^(eta-1))/(1-theta_1)) )^(1/(eta-1));
Psi_51   = ((1-theta_5*betta*pibar^(eta-1))/(1-theta_5*betta*pibar^(eta))) / ((1-theta_1*betta*pibar^(eta-1))/(1-theta_1*betta*pibar^(eta)));
if Omega_21<1e-6 || Psi_21<1e-6; error('something wrong!'); end
aux1 = ((1-theta_1*pibar^(eta-1))/(1-theta_1))^(1/(eta-1)) * eta/(eta-1)*(1-theta_1*betta*pibar^(eta-1))/(1-theta_1*betta*pibar^(eta));
aux2 = (1/5+1/5*(Omega_21*Psi_21)^(1-eta)+1/5*(Omega_31*Psi_31)^(1-eta)+1/5*(Omega_41*Psi_41)^(1-eta)+1/5*(Omega_51*Psi_51)^(1-eta))^(1/(eta-1));
mctildebar = aux2/aux1;

wbar = mctildebar;      % Real wage

% Reset prices
pstarbar_1 = eta/(eta-1)*(1-theta_1*betta*pibar^(eta-1))/(1-theta_1*betta*pibar^(eta))*mctildebar;
pstarbar_2 = eta/(eta-1)*(1-theta_2*betta*pibar^(eta-1))/(1-theta_2*betta*pibar^(eta))*mctildebar;
pstarbar_3 = eta/(eta-1)*(1-theta_3*betta*pibar^(eta-1))/(1-theta_3*betta*pibar^(eta))*mctildebar;
pstarbar_4 = eta/(eta-1)*(1-theta_4*betta*pibar^(eta-1))/(1-theta_4*betta*pibar^(eta))*mctildebar;
pstarbar_5 = eta/(eta-1)*(1-theta_5*betta*pibar^(eta-1))/(1-theta_5*betta*pibar^(eta))*mctildebar;

% Average prices
pbar_1     = ((1-theta_1*pibar^(eta-1))/(1-theta_1))^(1/(eta-1))*pstarbar_1;
pbar_2     = ((1-theta_2*pibar^(eta-1))/(1-theta_2))^(1/(eta-1))*pstarbar_2;
pbar_3     = ((1-theta_3*pibar^(eta-1))/(1-theta_3))^(1/(eta-1))*pstarbar_3;
pbar_4     = ((1-theta_4*pibar^(eta-1))/(1-theta_4))^(1/(eta-1))*pstarbar_4;
pbar_5     = ((1-theta_5*pibar^(eta-1))/(1-theta_5))^(1/(eta-1))*pstarbar_5;
test = 1/5*pbar_1^(1-eta) + 1/5*pbar_2^(1-eta) + 1/5*pbar_3^(1-eta) + 1/5*pbar_4^(1-eta) + 1/5*pbar_5^(1-eta) - 1;
if max(abs(test))>1e-6; error('something wrong!'); end

% Price dispersion
sbar_1 = (1-theta_1)/(1-theta_1*pibar^(eta))*pstarbar_1^(-eta);
sbar_2 = (1-theta_2)/(1-theta_2*pibar^(eta))*pstarbar_2^(-eta);
sbar_3 = (1-theta_3)/(1-theta_3*pibar^(eta))*pstarbar_3^(-eta);
sbar_4 = (1-theta_4)/(1-theta_4*pibar^(eta))*pstarbar_4^(-eta);
sbar_5 = (1-theta_5)/(1-theta_5*pibar^(eta))*pstarbar_5^(-eta);
sbar   = 1/5*sbar_1 + 1/5*sbar_2 + 1/5*sbar_3 + 1/5*sbar_4 + 1/5*sbar_5;

% Labor, output, consumption
nbar     = (wbar*sbar^siggma)^(1/(siggma+varphi));    
ybar     = nbar/sbar;
ybar_sss = ybar;
cbar     = ybar;

% Aggregate TFP, natural output 
tfpbar   = abar/sbar;
ynat     = ybar;
ynat_mis = ybar;

% Markups
markupbar_1 = pbar_1*mctildebar^(-1);
markupbar_2 = pbar_2*mctildebar^(-1);
markupbar_3 = pbar_3*mctildebar^(-1);
markupbar_4 = pbar_4*mctildebar^(-1);
markupbar_5 = pbar_5*mctildebar^(-1);

% Numerator and denominator of optimal reset price (G.2)
psibar_1 = mctildebar*ybar^(1-siggma)/(1-theta_1*betta*pibar^eta);
psibar_2 = mctildebar*ybar^(1-siggma)/(1-theta_2*betta*pibar^eta);
psibar_3 = mctildebar*ybar^(1-siggma)/(1-theta_3*betta*pibar^eta);
psibar_4 = mctildebar*ybar^(1-siggma)/(1-theta_4*betta*pibar^eta);
psibar_5 = mctildebar*ybar^(1-siggma)/(1-theta_5*betta*pibar^eta);
phibar_1 = ybar^(1-siggma)/(1-theta_1*betta*pibar^(eta-1));
phibar_2 = ybar^(1-siggma)/(1-theta_2*betta*pibar^(eta-1));
phibar_3 = ybar^(1-siggma)/(1-theta_3*betta*pibar^(eta-1));
phibar_4 = ybar^(1-siggma)/(1-theta_4*betta*pibar^(eta-1));
phibar_5 = ybar^(1-siggma)/(1-theta_5*betta*pibar^(eta-1));





